Increasing the density melts ultrasoft colloidal glasses 
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We use theory and simulations to investigate the existence of amorphous glassy states in ultrasoft 
colloids. We combine the hyper-netted chain approximation with mode-coupling theory to study the 
dynamic phase diagram of soft repulsive spheres interacting with a Hertzian potential, focusing on 
low temperatures and large densities. At constant temperature, we find that an amorphous glassy 
state is entered upon compression, as in colloidal hard spheres, but the glass unexpectedly melts 
when density increases further. We attribute this re-entrant fluid-glass transition to particle softness, 
and correlate this behaviour to previously reported anomalies in soft systems, thus emphasizing its 
generality. The predicted fiuid-glass-fluid sequence is confirmed numerically. 
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The hard sphere system is an emblematic model in sta- 
tistical physics. It is a useful theoretical model for liquids 
which can be studied experimentally using colloids 
its density driven fluid-crystal and fluid-glass transitions 
have been extensively studied. Recently, there has been 
growing interest in systems composed of ultrasoft parti- 
cles, which can be realized experimentally using colloidal 
particles made of star-shaped polymers or microgels Q. 
Since softness introduces (at least) one more dimension 
to the phase diagram, new physical behaviour can be 
expected from the competition between entropy, which 
governs the behaviour of hard spheres, and energy, which 
appears when soft particles overlap. 

As a result, the physics of systems consisting of soft 
particles is considerably richer than that of hard ones. 
Unusual equilibrium phase diagrams arise depending on 
the particular features of particle softness. These include 
cluster crystal structures at arbitrary large densities Q 
and a non-monotonic density dependence of the freez- 
ing transition with a complex cascade of crystalline 
states In the latter case, the thermodynamically sta- 
ble fluid phase also exhibits density anomalies, at the 
structural [sHU and dynamical @, H, E3] levels. 

Exploring the existence of amorphous solid states at 
low temperature and large density in soft particle mate- 
rials is an emerging field. The glassy dynamics of star 
polymers have been investigated both theoretically and 
experimentally, and, in particular, non-trivial phase di- 
agrams and novel glassy states have been described in 
star-polymer mixtures [2|, |llj, |l2|. Several studies were 



dedicated to the perturbations brought about by adding 
an energy scale to the hard sphere glassy behaviour. 
Such a perturbation not only shifts the glass transition 
to larger densities 13, UH but also affects qualitatively 
the dynamics near the transition, in particular the glass 
fragility IU IH 3]. Experiments have also investigated 



the physics of glass aging [17[, mostly in microgels, with a 
particular interest on the jamming transition intervening 
deep into the glass phase [3, [l9[ . 

In this work we investigate whether the anomalous 
equilibrium phase diagram and transport properties of 
soft repulsive sphere systems also impact their transition 
to a non-ergodic glassy state. This possibility was evoked 
for star polymers [llj for a narrow range of star function- 
alities, but was not confirmed numerically because crys- 
tallization intervened before slow dynamics set in. Other 
theoretical studies of soft repulsive spheres focused on 
densities close to the hard sphere transition, leaving the 
behaviour at large densities unexplored fLU ] . 

Our main finding is that an 'anomalous' density depen- 
dence of the glass transition temperature is found both 
theoretically and in computer simulations for the well- 
studied Hertzian model for ultrasoft particles. Thus, it 
is possible to melt a glass made of soft particles by in- 
creasing its density, a scenario that has not been observed 
yet in experiments. This is surprising on general grounds, 
because particle crowding typically yields jammed, rather 
than fluid, states for 'usual' particles such as hard spheres 
or molecular fluids. Indeed, traffic jams are rarely eased 
by putting more cars on the road. 

We consider N particles in a volume V interacting with 
a soft, repulsive, pairwise additive Hertzian potential, 



V(r < a) = e(l - rja) 



5/2 



(1) 



and V(r > a) = 0. Originally derived to describe 
the elastic repulsion between deformablc spheres, the 
Hertzian potential has recently been used to study theo- 
retically the physics of soft repulsive particles in various 
contexts, and seems to describe well the interaction be- 
tween microgel colloids jlq . Its physics is very similar 



to the harmonic sphere model inspired by thephysics of 
foams (joj . or to the Gaussian core model [ij, derived 



2 



from polymer physics. The control parameters of the 
system are the particle number density, p = Na 3 /V, and 
the ratio of the temperature T and energy scale e. Since 
particles can overlap at a finite energy cost, p can sig- 
nificantly exceed 1 and, indeed, we are interested in this 
high density regime. We use reduced units, measuring 
lengths in units of a, and temperatures in units of e. 

We study the structure and dynamics of Hertzian 
spheres combining theory and simulations. For the the- 
ory, we closely follow an earlier investigation of har- 
monic spheres to which we refer for further details and 
equations [l4| . We examine a monodisperse system of 
Hertzian spheres using the hyper-netted chain (HNC) ap- 
proximation 2l| to obtain the pair correlation function 
g(r) and the structure factor S(q) for a given state point 
(p,T). We then inject S(q) in the dynamical equations 
of motion for time correlation functions derived within 
mode-coupling theory (MCT) [22j , which we solve numer- 
ically [H]. We obtain the time dependence of density- 
density correlation functions for a wide range of wave- 
vectors for a broad range of densities and temperatures. 

To verify theoretical predictions we perform molecular 
dynamics (MD) simulations of a polydispersc Hertzian 
sphere system. Polydispersity is needed to prevent 
crystallization or fractionation that easily occur in the 
monodisperse system. The particle diameters Oi are ran- 
domly drawn from a flat probability distribution with 
0.826 < Gi < 1.126. Two particles i,j interact through 
Eq. (P with a = {a, + aj)/2. We use N = 9000 and 
periodic boundary conditions. Equilibration is achieved 
through periodic velocity rescaling, production runs are 
carried out in the microcanonical ensemble. 

We now discuss two representative quantities charac- 
terizing the structure of the fluid at various state points. 
In Fig. [T}a, we show the evolution with the density of the 
pair correlation function g(r) at constant temperature, 
T = 0.00287, as predicted by the HNC approximation. 
As expected, with increasing density particles are getting 
closer to each other and as a result, the position of the 
first peak of g(r) moves towards smaller values of r. The 
monotonic shift of the peak position of <?(r) is accompa- 
nied by an 'anomalous', non-monotonic evolution of the 
peak height Q ■ Specifically, with increasing density from 
p = 1.1 to 1.7 the height of the first peak increases, but it 
decreases when increasing density further from p = 1.7 to 
2.84. This evolution is 'anomalous' in the sense that the 
first peak of g{r) typically increases monotonically with 
density in simple liquids at constant temperature. The 
locus of the density anomaly evolves with temperature, 
as studied in detail for the harmonic potential |8( . 

Since g(r) and the structure factor S(q) are related via 
Fourier transform (2lj , it is not surprising that the same 
anomalous trend is observed in Fig. [T}c, which shows 
the HNC-derived S(q) for the same temperature. With 
increasing density the position of the first peak shifts 
to larger g, corresponding to smaller interparticlc dis- 
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FIG. 1: Nonmonotonic evolution of structural quantities with 
increasing density (arrows) as predicted by theory (top) and 
obtained from MD simulations (bottom), a) Predicted pair 
correlation function for T = 0.00287 and p = 1.1, 1.4, 1.7, 2.1, 
2.43, 2.84; b) Numerically measured g(r) for T = 0.0043 and 
p = 1.2, 1.45, 2.05, 3.0, 4.05; c) Predicted structure factor for 
T = 0.00287 and p = 1.1, 1.5, 2.1, 2.63, 2.84; d) Numerically 
measured S(q) for T = 0.0043 and p = 1.2, 1.7, 2.4, 3.0, 4.05. 



tances. Simultaneously, the peak height exhibits a non- 
monotonic variation. The maxima of g(r) and S(q) typ- 
ically occur at slightly different densities, respectively 
p « 1.7 and p « 2.6 at T = 0.00287. 

The trends exhibited by g(r) and S(q) calculated using 
the HNC approximation are also evident in the simula- 
tions. In Fig.[T]-b,d we show the pair correlation function 
and the structure factor obtained in simulations at con- 
stant T = 0.0043 and various densities. Both quantities 
exhibit a non-monotonic evolution with density, g(r) hav- 
ing its maximum near p k, 2.05, and S(q) near p « 2.4. 
Thus, the two maxima are closer to each other in sim- 
ulations than according to theory. However the simi- 
larity between theory and simulations is quite striking 
from these plots, showing that HNC is performing rea- 
sonably well for the explored state points. Quantitative 
agreement cannot be expected as we are using different 
particle size distributions in both approaches. In partic- 
ular, we believe that polydispersity is responsible for the 
fact that the first peaks of both g{r) and S{q) obtained 
from MD are quite a bit smaller and broader than those 
calculated using the HNC approximation. 

The non-monotonic dependence of the height of the 
first peak of g(r) on density was previously discussed for 
several similar models. Physically, it originates from a 
competition between energy and entropy with an 'anoma- 
lous' outcome due to the extreme softness of the parti- 
cles d-Q. Briefly, at smaller densities, the system can 
easily minimize its free energy by avoiding the energet- 



ically costly particle overlaps. In this regime, the sys- 
tem effectively behaves as hard spheres, with g(r) and 
S(q) becoming sharper as p increases. At larger densi- 
ties, it becomes difficult, and thus entropically costly, for 
the system to find configurations with no overlap (note 
that the nominal volume fraction <fi = ■np/& may exceed 
100%). Rather than trying to do the impossible, parti- 
cles increase their overlaps, have effectively more space, 
and may thus access a larger number of configurations. 
This increasing disorder leads to the broadening of the 
first peak of correlation functions. The entropy gain can 
be larger than the energy cost of the overlaps because 
the Hertzian potential is bounded at any separation. 

The above described competition between entropic and 
energetics effects is purely statistical. It is reasonable to 
expect that the resulting structural anomalies will cause a 
non-trivial density dependence of the equilibrium freezing 
transition Q. However, this competition and the anoma- 
lies are not obviously relevant to the nonequilibrium glass 
transition, which is related to the increasing difficulty to 
transport matter with increasing particle crowding. Our 
main result is that the anomalous density dependence of 
equilibrium structural quantities is indeed accompanied 
by the non-monotonic dependence of structural relax- 
ation and ergodicity breaking. Anomalous density de- 
pendence of transport properties were found before in 
the thermodynamically stable, non-glassy soft particle 
fluid 0, 0, [SHlij], but remained essentially unexplored 
for the deeply supercooled state and for the glass. 

In Fig. [2] we show the p dependence of the density- 
density time correlators as obtained from both theory 
and simulations. As for structural quantities we show the 
correlators at constant T and increasing p, which exposes 
the most remarkable features of the data, but a much 
broader range of state points was explored. Figures [2] 
exhibit the same qualitative behavior. With increasing 
density the relaxation becomes progressively slower and 
the fluid approaches a glass transition qualitatively sim- 
ilar to the one observed with hard particles. A two-step 
decay of time correlation functions emerges, with the slow 
relaxation being strongly dependent on density. For the 
temperatures shown in Fig. [2j F(q,t) does not decay to 
zero above p re 1.9 within MCT at T = 0.00287, and 
above p re 2 within MD at T = 0.003. Although the 
relaxation time is infinite within theory, we can simply 
report that it has become too large to be measured in the 
simulations and that, for practical purposes, the system 
is a non-crgodic glass. In the glass phase there is, at first, 
little change when p increases, apart from a weak varia- 
tion of the long-time limit of the correlators. However, 
quite surprisingly, the glass suddenly melts when it is 
compressed further, and the relaxation time becomes fi- 
nite again. This is counter-intuitive from the standpoint 
of glass physics, although perhaps not totally unexpected 
for a system known to displaying other anomalies. Past 
the glass melting transition, the dynamics becomes faster 




FIG. 2: Nonmonotonic evolution of time correlation functions 
with density at constant temperature for theory (top, q = 
9.1) and simulations (bottom, q = 9.4). A two-step decay 
initially develops when p increases, non-ergodic glass states 
are obtained at intermediate densities, and the glass melts 
upon further compression. 



when density is increased even further. 

The above description of time-dependent correlators 
suggests the existence of a reentrant glass transition line 
in the (p, T) phase diagram, which we confirm in Fig. [3] 
In the theory we easily obtain the critical mode-coupling 
temperature, T c (p), with arbitrary precision, by fitting 
the temperature evolution at constant density of the 
structural relaxation time to a power law divergence, 
r oc (T — T c ) -7 [22|]. A finite mode-coupling temper- 
ature emerges in a way qualitatively similar to the re- 
sults obtained for the harmonic potential |14j . However, 
the mode-coupling transition line has a clear maximum, 
T c re 0.003, near p re 2.4, slightly below the density where 
S(q) has its maximum. It then decreases smoothly with 
increasing p. When performed for simulation data, the 
same analysis works well over a moderate range of 2- 
3 decades of r but, as usual [22j| . the MCT power law 
breaks down below some T. Remarkably, however, the 
numerical T c values reported in Fig. [3j exhibit trends in 
good qualitative agreement with theory. The important 
point is the robustness of the reentrant glass transition 
line with a maximum around p re 2.4, at nearly the same 
density where the freezing transition temperature peaks 
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FIG. 3: Dynamic phase diagram showing the exact (theory) 
and fitted (simulations) critical temperatures obtained from 
power law fits to the structural relaxation times. A reentrant 
glass transition is obtained in both cases with a maximum 
near p ~ 2.4. 



for the monodisperse system Q, T m ps 0.008. 

The quantitative disagreement between the two glass 
lines has multiple origins. First, the theoretical results 
were obtained for the monodisperse system, but the simu- 
lations were performed using a polydisperse one. Second, 
both HNC and MCT have their well known limitations. 
It is possible to improve the static approximation plj . 
but the predicted (spurious) divergence of the relaxation 
time predicted by MCT will remain. Moreover, MCT 
predictions stem from the knowledge of pair correlation 
functions, incorrectly ignoring additional structural in- 
formation [24J]. Nevertheless, in the present case, the 
evolution of pair correlations are certainly the major ori- 
gin of the reentrant glass transition, which is thus suc- 
cessfully captured by MCT. Analogous successes of MCT 
occurred in the past, for instance for sticky colloidal par- 
ticles [25[. Despite its limitations, our procedure should 
therefore be useful in revealing new dynamical features 
relevant in the context of colloidal physics with soft in- 
teractions, as our numerical simulations confirm. 

In summary, by combining simple analytical ap- 
proaches and simulations, we predicted and observed nu- 
merically the existence of a reentrant fluid-glass transi- 
tion in a system of Hertzian spheres at large densities. We 
believe that our results are not restricted to this particu- 
lar interaction and similar behavior should arise in other 
systems with soft repulsive potentials, and should find 
experimental realizations using ultrasoft colloids. Our re- 
sults suggest that making softer colloids not only affects 
the kinetic glass fragility H, 15, 3], but can also pro- 



duce new, non-intuitive effects such as a reentrant glass 
transition, which remains to be seen in an experiment. 
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